Dipolar Bose-Einstein condensates with dipole dependent scattering length 
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We consider a Bose-Einstein condensate of polar molecules in a harmonic trap, where the effective 
dipole may be tuned by an external field. We demonstrate that taking into account the dependence 
of the scattering length on the dipole moment is essential to reproducing the correct energies and for 
predicting the stability of the condensate. We do this by comparing Gross-Pitaevskii calculations 
with diffusion Monte Carlo calculations. We find very good agreement between the results obtained 
by these two approaches once the dipole dependence of the scattering length is taken into account. 
We also examine the behavior of the condensate in non-isotropic traps. 



INTRODUCTION 



Degenerate atomic quantum gases are typically very 
dilute systems. Nevertheless, inter-atomic interactions 
strongly determine many of the observed phenomena and 
their underlying physics Q, Q ■ Until recently, only short- 
range and isotropic interactions have been considered. 
However, recent developments in the manipulation of 
cold atoms and molecules have been paving the way to- 
wards the analyses of polar gases in which dipole-dipole 
inter-particle interactions are important 0, 0>M IE 01 0- 
A major break through has been very recently achieved 
by the experimental group in Stuttgart 0, where a 
Bose-Einstein condensate (BEC) of strongly magnetic 
52 Cr atoms has been realized. This experiment ob- 
served the effects of dipole-dipole interactions on the 
shape of the condensate. New exciting phenomena are 
expected to occur in these quantum gases since the par- 
ticles interact via dipole-dipole interactions which are 
long-ranged and anisotropic. Recent theoretical analy- 
ses have shown that the stability and excitations of dipo- 
lar gase s are crucially determined by the trap geometry 
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Yi and You [Til IT^ | first introduced a pseudopotential 
appropriate for describing slowly moving particles inter- 
acting via short range repulsive forces and long range 
dipolar forces. The dipolar long range part of their pseu- 
dopotential is identical to the long range part of the orig- 
inal potential. The pseudopotential also includes a con- 
tact (delta function) potential whose coefficient is propor- 
tional to the scattering length. For non-polar particles, 
the scattering length is solely due to short range inter- 
actions. A crucial point that Yi and You have shown 
is that the scattering length is also dependent on the 
long-range dipolar interaction. They have argued for the 
need to take this into account when calculating conden- 
sate properties through the Gross-Pitaevaskii equation 
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(GPE). However, they stopped short at actually explic- 
itly showing how condensate properties are influenced by 
the dependence of the scattering length on the dipole mo- 
ment. Rather, they have analyzed the dependence of the 
condensate properties on the ratio of the dipole moment 
to the scattering length. This is appropriate if one con- 
siders a fixed dipole moment (and scattering length) but 
varies the number of particles in the trap. On the other 
hand, the scenario we wish to consider is that where dipo- 
lar interactions are tuned by an external field, in which 
case the dipole moment and scattering length do not scale 
equally. Also, Ref. [12| did not consider negative dipolc- 
dependent scattering lengths. 

In this work we focus on a trapped gas of dipolar par- 
ticles, where the inter-particle interaction is dominated 
by the dipole-dipole force. A possible realization includes 
an (electrically polarized) gas of heteronuclear molecules 
with a large permanent electric dipole moment. The ef- 
fective dipolar interaction may be tuned by the competi- 
tion between an orienting electric field and the quantum 
or thermal rotation of the molecule. For example, the 
2 n 3 / 2 ground state of the OH molecule is completely po- 
larized in a field of about 10 4 V/cm. For smaller fields, 
the field strength determines the degree of polarization 
and the size of the dipole moment. Another approach 
for tuning the dipolar interaction, by using a rotating 
external field, was proposed in Ref. |19|. 

We consider the case of a bare (zero dipole mo- 
ment) positive scattering length, and show that as the 
dipole moment is increased, this scattering length be- 
comes smaller and, above a certain dipole strength, neg- 
ative, and then again positive after crossing a resonance. 
Taking this variation of the scattering length into ac- 
count is necessary to describe the creation of new two- 
body bound states. This has consequences for the ba- 
sic theory, and for the predicted stability of BECs with 
anisotropic interactions, especially for polar molecules, 
which can have large dipole moments. To verify that the 
dipole-dependent scattering length is essential to repro- 
ducing the correct energetics, and to ascertain the va- 
lidity of this approach, we have solved the many-body 
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Schrodinger equation using the diffusion Monte Carlo 
(DMC) method, as well as solved the Schrodinger equa- 
tion by direct diagonalization for two particles. In these 
calculations we used a potential consisting of short range 
hard wall potential and long range dipolar interaction. 
We then compare these results with solutions of the GPE, 
in which we employ the appropriate pseudopotential with 
the dipole-dependent scattering length. 

Section [H] describes the system under study while 
Sec. lllll presents our results for dipolar gases under spher- 
ically symmetric and cylindrically symmetric confine- 
ment. Selected results for systems under spherically sym- 
metric traps have been presented in a previous paper [2(| . 
Here, we provide numerical details (see the Appendix), 
justify the formalism employed in detail, and significantly 



extend the discussion of our findings. Section IIVI con- 
cludes. 



II. FORMULATION 



For N identical bosons in an external trap potential 
Vext{ r ) with pair-wise interaction V(r — r'), at zero 
temperature, the condensate is typically described us- 
ing mean-field theory. All the particles in the condensate 
then have the same wave function ip(r). As a first at- 
tempt, one may write a Hartree-Fock equation for this 
wavefunction |2lj : 



J 



-^ 2 + V ext {r) + (N 
2m 



1) / dr'V{r-r')\tp(r')\ 



ip(r) 



(1) 



where fj, denotes the chemical potential, r the displace- 
ment from the trap center, and m the mass of a particle. 
In the following we restrict ourselves to cylindrically sym- 
metric harmonic traps with V ext {r) = ^m{L0 2 p p 2 + uj 2 z 2 ). 

We consider an interaction potential V(r) which con- 
sists of a dipolar interaction and a short range hard wall 
with cutoff radius b: 



V(r) = 



d 

oc 



2 1-3 cos^ 9 



if r > b 
if r < b 



(2) 



where d denotes the dipole moment (in Gaussian units), 
r the distance vector between the dipoles, and 9 the angle 
between the vector r and the dipole axis, which we take 
to be aligned along the z-axis of the trap. The hard wall 
cutoff corresponds to hard particles with diameter b (two 
particles cannot penetrate each other when the distance 
between their centers is equal to their diameter). Unfor- 
tunately, the short range part of the interaction potential 
V causes the integral in Eq. Q to diverge. Fortunately, 
this divergence can be cured since the condensate is very 
dilute and ultra-cold, which implies that the particles in 
it are moving very slowly. Thus the potential V may be 
replaced by an effective potential (pseudopotential) V e f / 
with a milder small r behavior, which reproduces the 
two-body scattering wavefunction asymptotically in the 
zero energy limit, and which does not lead to divergencies 
when used in the Hartree-Fock mean-field description. 

For a short range potential V, the low energy scattering 
amplitude is completely determined by one parameter, 
the scattering length a, and an appropriate pseudopo- 
tential is V e ff(r) — 4 ^ a S(r). In general, the pseudopo- 
tential is chosen such that its first-order Born scattering 
amplitude reproduces the complete scattering amplitude 
of the original potential V in the zero-energy limit. For 



a potential with long range dipolar part, Yi and You 
[Iull2| proposed the following pseudopotential: 



AirTi 2 a(d) 
Veff(r) = S(r) 



.1-3 cos 2 e 



(3) 



where the scattering length a(d) depends on the dipole 
moment d. Note that the long range part of the pseu- 
dopotential is identical to the long range part of the orig- 
inal potential. By construction, the scattering amplitude 
of V e f f calculated in the first-order Born approximation 
agrees with the full zero-energy scattering amplitude of 
V. To verify the validity of the pseudopotential for sys- 
tems of experimental interest, we compute the low energy 
scattering amplitude f(k,k') for two OH-like molecules 
interacting through the model potential given in Eq. J2J), 
and compare it with the scattering amplitude computed 
in the first Born approximation. 

The low energy scattering amplitude in the presence 
of the long-range dipolar interaction may be expanded in 
partial waves: 

f(k,k')=in J2 *tm' {k)YC m {k)Y Vm ,{k'), (4) 

IrnA'm' 



with t\™ (fc) the reduced T-matrix elements. The 
anisotropic dipolar potential implies that / depends on 
the incident and scattered directions k and k' . The 
reduced T-matrix elements are related to the usual T- 
matrix elements T}^ 1 ' (k) = (lm\T(k)\l'm') by tf™'(k) = 



(fc) 



— 2fe — ■ F° r fc — > they are energy independent, and 
act as generalized scattering lengths. In particular, the 
scattering length a(d) is given by too(0) and depends on 
the dipole moment d. 
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FIG. 1: Absolute values of selected reduced T-matrix ele- 
ments tni = t\o (symbols) for the OH-OH model potential as 
a function of the relative scattering energy, compared with the 
first-order Born approximation (dashed lines). Note, the first- 
order Born approximation to too diverges and is not shown. 



Figure ^ compares some reduced T-matrix elements 
(symbols), computed through numerical close-coupling 
calculations, with the first Born approximation [ill 1 1 21 
(dashed lines) , for the potential given in Eq. (0 with 
parameters chosen to describe the scattering between two 



rigid OH-like molecules. In particular, we have taken 
to = 17.00 amu, d = 0.66 a.u., and b = 105 a.u. (this is a 
reasonable value for a molecular scattering length). For 
the <qq channel, the Born approximation to the long range 
dipolar part of the potential gives no contribution, while 
that to the hard-core part diverges (and is therefore not 
shown in the figure). For the other channels, there is a re- 
markably good agreement in the E — > limit between the 
exact reduced zero-energy T-matrix elements and those 
calculated in the first Born approximation. The agree- 
ment becomes less good at finite but small scattering en- 
ergies (of the order 10 -7 K). Also, too(k) is not constant 
(as it would be in the threshold limit) even for very low 
energies of the order of 10~ 10 K. This suggests that, e.g., 
an effective range correction pli. l24j may become impor- 
tant at finite E, but here we consider only the E = 
limit. The £qq (k = 0) matrix element determines a(d), 
the dipole-dependent scattering length of the pseudopo- 
tential, Eq. lyj. In this way, the Born approximation to 
the pseudopotential gives the correct ^(k = 0) value. 
The fact that a(d) depends on d is important, since the 
strength of the dipolar interactions may be controlled by 
an external field. Our analysis confirms that the pseu- 
dopotential approximation provides a good description 
in regions away from resonance (see also below) . 

Replacing V in Eq. Q with V e ff of Eq. with 
the scattering length a(d) determined through numerical 
coupled-channel calculations for the model potential V, 
we obtain the (time-independent) Gross-Pitaevskii equa- 
tion (GPE): 



fiijj(r) 



-|^V 2 + V ext (r) + (N-l)[ dr'V ef f(r - r')|^(r')| 2 



V>(r). 



(5) 



For the following, we define a dipole length = 
md 2 /h 2 . This is the distance at which the dipolar po- 
tential energy equals the kinetic energy (estimated from 
the uncertainty relation) of two interacting dipoles. In 
Fig. [2] we plot the ratio a/b of the scattering length a to 
the hardwall cutoff b as function of D*/b. This provides 
a universal curve for the model potential given in Eq. J2J) , 
which determines the scattering length for any given cut- 
off b and dipole length Z?» . Note the appearance of res- 
onances, corresponding to the appearance of new bound 
states. In the neighborhood of a resonance, the scatter- 
ing length tends towards — oo before and +00 after the 
resonance. These resonances have been identified before 
in dipolar scattering 0, [2f| . They would occur at fields 
of order MV/cm in atoms and kV/cm in heteronuclear 
molecules, or perhaps at 10 5 V/cm in atoms, if assisted 
by Feshbach resonances [26|. 

It is instructive to connect Fig. [21 back to concrete 
dipoles that can be handled in the laboratory. Consider, 



for example, atomic chromium, which is generating a lot 
of interest now that it has been Bose condensed. 52 Cr has 
a magnetic dipole moment of 6/13, and a sextet scatter- 
ing length of 112 a.u. Identifying the hardwall cutoff with 
this scattering length, chromium would appear on Fig. [21 
at the value D*/b = 0.4, and the scattering length would 
be renormalized by ~ 0.6% of its value in the absence 
of a dipole moment. Indeed, this correction is already 
included in two-body modeling of the Cr-Cr interaction 
[271 |28| . For chromium, resonances of the type shown in 
Fig. [21 play no role. 

For a heteronuclear polar molecule, however, the situ- 
ation can be quite different. Consider, for example, the 
OH radical, which is also the focus of intense experimen- 
tal efforts 0, |2{j . This molecule has a permanent electric 
dipole moment of 0.66 atomic units, and therefore a huge 
dipole length, = 1.35 x 10 4 a.u.. If we assume a small 
cutoff, such as b — 105 a.u., then D^/b = 128, which 
is way off-scale in Fig. [21 In other words, this potential 



4 




10 15 

DJb 



LU 

w 
"E 



LU 



1.6 
1.5 
1.4 
1.3 
1.2 



1.1 









+ Monte-Carlo 

Gross-Pitaevskii 

Constant scattering length 

Contact potential only 









0.04 



0.08 



Dipole length (units of a ho ) 



0.12 



FIG. 2: Scattering length a(d) versus dipole length D, for the 
dipolar potential with hard wall cutoff b given in Eq. (|5J . 



supports a large number of OH-OH dimer states. These 
dimer states could play a large role in the physics of a 
gas composed of polar molecules. 



FIG. 3: Energy per particle E/N as a function of the dipole 
length D, for N = 10. The symbols with error bars show 
DMC results. The solid line is the solution of the GPE. The 
dotted line shows a GP calculation without taking the dipole 
dependence of the scattering length into account, setting it 
equal to the hardwall cutoff. The dashed line shows a GP 
calculation with the dipole dependent scattering length, but 
omitting the long range dipolar term. The lines terminate at 
the collapse point of the condensate. 



III. MANY BODY AND TWO-BODY 
CALCULATIONS 

The mean field GP approach to solving the many-body 
dynamics is necessarily approximate. However, the full 
iV-body Schrodinger equation with the actual interaction 
potential V, Eq. can be solved numerically for rela- 
tively small N. For N — 2, the Schrodinger equation can 
be solved by direct diagonalization. For N > 2, direct 
diagonalization becomes impractical and we solve the 
Schrodinger equation instead by the DMC method [3u| . 
This section compares the results obtained by the three 
methods. Details of how these schemes are implemented 
are presented in the Appendix. 

We note that the system is scalable: if the hardwall 
cutoff b and the dipole length are each scaled by a 
factor K, then the scattering length scales by the same 

amount. If the harmonic oscillator lengths \/—^— and 

° y mu z 

\Jj£j- are also scaled accordingly, then the entire spec- 
trum (both the exact spectrum obtained by solving the 
iV-body Schrodinger equation and the spectrum obtained 
within the mean field approximation) remains the same 
apart from scaling by a factor 1/K 2 . 

In what follows, unless stated otherwise, we work with 
an isotropic harmonic trap: u> p — u> z = ui. Correspond- 
ingly, the natural unit of length is the harmonic oscilla- 
tor length aho = *J ^j- In the following, we consider the 
two-body potential V, Eq. @, for two different values of 
b, i.e., b — 0.0137 a.u. and b — 0.0433 a.u., and varying 
d. For concreteness, OH molecules in a trap with a fre- 



quency uj — 2n kHz would have a^o = 5800 a.u. so that 
b = 0.0137a/! O corresponds to a hardwall cutoff of 79 a.u. 



A. Energies and collapse in an isotropic trap 

In this paper, our primary interest is in BEC's of polar 
molecules. In such systems, the dipole moment is some- 
thing that is directly under the experimentalists control: 
in zero electric field the dipole moment vanishes, while at 
higher fields its value can be continuously tuned. For this 
reason, we consider properties of dipolar condensates as 
a function of the dipole moment, which is taken as a sub- 
stitute for the dependence on the strength of an external 
electric field. One quantity can be converted into the 
other, of course, via the polarizability of the molecule. 

Figure [3] presents the ground state energy per parti- 
cle in units of the harmonic oscillator energy Eho = Tv+> 
versus the dipole length in units of the harmonic oscil- 
lator length aho for N = 10 molecules. In this figure, the 
solid line shows the GP energies obtained using V e f / with 
the dipole-dependent s-wave scattering length a(d). For 
comparison, the symbols show the results from the DMC 
simulations, which solve the iV-body Schrodinger equa- 
tion for the model potential V, Eq. J2J (statistical uncer- 
tainties are indicated by vertical errorbars). The agree- 
ment is excellent, attesting to the validity of parametriz- 
ing the GP equation in terms of the dipole-dependent 
scattering length a(d). 

For zero dipole moment, the scattering length is equal 
to the hardwall cutoff b. Thus the energy per particle 



5 



is larger than its noninteracting value of 1.5K.W. As the 
dipole moment increases, the energy per particle drops. 
Qualitatively, this has been explained previously by an 
elongation of the dipolar gas in the z-direction. This 
allows more dipoles to encounter one another in an at- 
tractive "head-to-tail" orientation. As shown below, we 
indeed observe an elongation of the condensate. 

However, by far the greatest influence of the dipolar 
interaction on the condensate energy comes from the 
dipole dependent scattering length. To illustrate this, 
the dashed line in Fig. [5] shows the result of a GP calcu- 
lation, which includes the contact (delta function) term 
oiV e ff with the dipole dependent scattering length, but 
omits the non-isotropic dipolar long range term. It is 
quite surprising how good this approximation is, in which 
the condensate is completely isotropic. For comparison, a 
dotted line shows the results from a GP calculation which 
keeps the long range interaction, but does not take the 
dependence of the scattering length on the dipole mo- 
ment into account. Instead, we use the zero-field scatter- 
ing length, i.e., a = b. It is clear that this approximation 
greatly overestimates the energy. It is not even a very 
good approximation for small dipole moments (the sec- 
ond derivative at d = does not match the curvature for 
the "correct" GP solution (solid line)). 

As the dipole moment increases, the condensate 
quickly moves toward a complete collapse, with E/N — » 
-co. This occurs because the scattering length attains a 
large, negative value. It is well known that a condensate 
with short range interaction only collapses when a critical 
combination of particle number and scattering length is 
reached, (JV-l)a = -0.57497a fco [sHIH. This is exactly 
the point where the dashed line in Fig. [31 terminates. It 
is also very nearly the point where the solid line, which 
provides the most complete description of dipolar gases 
at the GP level, terminates. The collapse discussed here 
is very different from that usually discussed in the con- 
text of dipolar gases, which takes the scattering length 
to be constant. In the latter case, collapse can also occur 
due to the presence of a large dipole moment, but (as 
seen from the dotted line in Fig. \5§ this collapse occurs 
for much larger dipole moments than the value predicted 
here by taking the dependence of the scattering length on 
the dipole moment into account. We therefore predict a 
kind of electric-field-induced condensate collapse, which 
shares many similarities with the collapse encountered 
in alkali atom BEC's near magnetic field Feshbach res- 
onances [82L Is^ . In particular, ramping the field across 
a resonance would likely result in the creation of dimer 
states of the molecules. 

Figure 0] shows the dependence of the energy per par- 
ticle on the dipole moment for various particle numbers, 
TV = 10, 20 and 50. In all cases the conclusions are the 
same. Namely, the GP energies obtained using V e // with 
the renormalized scattering length (solid lines) are an ex- 
cellent approximation to the DMC energies (symbols). 
Also, the collapse occurs for smaller dipole moment as 
the particle number increases, in accord with the collapse 
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FIG. 4: Energy per particle E/N as a function of the dipole 
length D, for TV = 10, 20 and 50 (from top to bottom). Solid 
lines show mean field GP energies. Dots (TV = 10), circles 
(TV = 20) and triangles (TV = 50) show DMC energies. Error 
bars indicate the statistical uncertainties of the DMC ener- 
gies. 
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FIG. 5: Variational many-body energies Ev/N for TV = 20, 
b = 0.0137a ho and D, = 0.079a ho (crosses), D, = 0.084a ho 
(squares) and D, — 0.093ah o (triangles) as a function of the 
Gaussian width b r , where b r — b p = b z . 



criterion. 

In the present discussion we have tuned the dipole mo- 
ment starting from d — 0. Because of the resonance 
structure shown in Fig. |3 it is possible to start with 
a large dipole moment and a small dipole-dependent 
scattering length. In this case we expect that the 
anisotropic long-range part of the effective potential will 
carry greater weight (i.e., will have a larger impact on 
the energetics and collapse). We have not considered 
this case in the present work. 

We now investigate the stability of dipolar Bose gases 
by considering the variational energy Ey obtained by a 
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variational many-body ansatz (see the Appendix for de- 
tails and notation) . Symbols in Fig. show Ey /N for 
b = 0.0137<2fj O , N — 20 and three different values of D*, 
i.e., = 0.079aho (crosses), £>» = 0.084afc o (squares) 
and D* = 0.0930^0 (triangles), as a function of the Gaus- 
sian width b r , which is treated as a variational param- 
eter (we set bp — b z — b r in this stability study). The 
variational parameters of F (Eq. I|A.6|> ) are optimized 
for — 0.079a/j O and 0.084a; lo by minimizing the en- 
ergy Ey for fixed b p and b z , i.e., for b p = b z = 1. For 
= 0.093a^ o , we use the same values for pi through p 5 
as for D* = 0.084a/j O since no local minimum exists for 
= 0.0930^0 (see below). 

Figurc[5]indicates that the variational energy for D* = 
0.079a/j O shows a local minimum at b r w la^o and a 
global minimum at b r w 0.09dh o . The "barrier" at 
b r w 0.3a; lo separates the large b r region in configura- 
tion space, where the metastable condensate state exists, 
from the small b r region where bound many-body states 
exist |42j. At very small b r , the variational energy be- 
comes large and positive due to the hardcore repulsion 
of the two-body potential. Figure [S] indicates that the 
energy barrier decreases with increasing . The dipolar 
gas collapses at the -D*/a/ lo value for which the energy 
barrier vanishes. Our DMC calculations show that the 
condensate prior to collapse is only slightly elongated, 
which justifies that our stability analysis parametrizes 
the one-body term 4> (see Eq. in terms of a single 

Gaussian width b r and which is consistent with our find- 
ing that the collapse is induced primarily by the negative 
value of a(d). We find similar results for other N values. 

We emphasize that the presence of the energy barrier 
is crucial for our DMC calculations to converge to the 
metastable condensate state for sufficiently large D^/aho 
and not to the cluster-like ground state (see also the Ap- 
pendix) . 



B. Condensate sizes and shapes 

The nature of the collapse is seen in more detail by 
looking at the size and shape of the condensate. To this 
end Fig. shows the root-mean-square widths Z and X 
of the condensate in the z and x directions, for A = 10 
and b = 0.01370/^ (the same parameters as in Fig. |3Jl. 
The solid and dashed lines are computed using the GP 
equation with the renormalized scattering length, and the 
symbols with error bars show, as before, the results of 
the DMC calculations; again, the agreement is excellent. 
For zero dipole moment, the condensate is isotropic and 
slightly larger than the harmonic oscillator length. As 
the dipole moment is turned on, the condensate contracts 
slightly in the x direction, and expands in the z direction. 
This illustrates the elongation that has been predicted for 
a dipolar condensate [l2fl . 

If in the effective potential we retain only the dipole- 
dependent scattering length, but omit the long range 
dipolar term, we obtain the dash-dotted line. This ap- 
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FIG. 6: Condensate sizes X = V < x 2 > and Z = V < z 2 >, 
for 10 particles with hard wall cutoff b = 0.0137ah o in a 
spherical trap. Solid and dashed lines show X and Z ob- 
tained by solving the GP equation using V e ; / with the dipole- 
normalized scattering length. Symbols with error bars show 
the corresponding DMC results. The dotted lines show results 
of a GP calculation with constant scattering length a = b. 
The dash-dotted line shows results of a GP calculation with 
dipole dependent scattering length that omits the long range 
dipolar interaction. Inset: The aspect ratio Z/X for the so- 
lution to the GP equation using V e ff (solid line), for the so- 
lution to the GP equation with constant scattering length 
a = b (dotted line), and for the DMC solution to the iV-body 
Schrodinger equation (symbols with error bars). 



proximation describes an isotropic condensate (Z = X) 
whose size is in-between Z and X obtained by the more 
complete description (dashed and solid lines). On the 
other hand, if we retain the dipolar long range term, but 
set the scattering length to a constant (its value for zero 
dipole moment, a — b), we obtain the dotted lines. This 
formulation overestimates the size of the condensate, but 
provides quite a good approximation to the aspect ratio 
Z/X (inset in Fig.^J - at least for small dipole moments 
As the dipole moment is increased, the condensate ap- 
proaches the point of collapse, and the aspect ratio cal- 
culated from the solution to the full GP equation (solid 
line in the inset) and the DMC calculations (symbols 
in the inset) increases rapidly. The approximation of a 
constant scattering length does not describe this behav- 
ior correctly. Interestingly, relatively large anisotropy is 
predicted just prior to collapse, even though the collapse 
mechanism is primarily s-wave dominated. 



C. Energetics for dipoles in non-isotropic traps 

We have also explored the behavior of dipolar BECs 
in a pancake-shaped trap with lu z /lu p = 10 and in a 
cigar-shaped trap with uj z /ujp = 1/10. For the pancake 
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FIG. 7: Energy per particle E/N obtained by solving the 
GP equation using V e ff as a function of the dipole length 
D* for N = 10 in a pancake-shaped trap with uj z /u] p = 10 
(dashed line) and in a cigar-shaped trap with lu z /u p =0.1 
(solid line). Symbols with error bars show the correponding 
DMC results, a^o is defined by the shortest dimension in each 
case (see text). The hard wall cutoff is b — 0.0137ah o . 



trap, we define = while for the cigar trap, 



aho = ^J^tr- he, compared to the spherical trap of 

the previous sections, we elongate the trap in either the 
axial or transverse directions. Figure [7\ indicates very 
good agreement between the energies calculated from 
the GP equation (solid and dashed lines) and from the 
TV-body Schrodinger equation (symbols). In the cigar- 
shaped trap, a non-vanishing dipole moment leads to a 
decrease of the energy compared to the energy obtained 
for d = 0. For the pancake-shaped trap, in contrast, the 
energy per particle increases for small dipole moments, 
and then decreases for larger dipole moments. In both 
traps, collapse occurs at the value at which the lines 
in Fig. [7| terminate. 

We note in particular the collapse of the condensate 
in the pancake trap. In [T^ | it was proposed that in 
a trap of such aspect ratio, if the scattering length is 
zero, there should never be collapse even for an arbitrar- 
ily large dipole moment. In our case, the "bare" (zero 
dipole moment) scattering length is positive, and so, ne- 
glecting the scattering length dependence on the dipole 
moment, the BEC should be even more resistant to col- 
lapse. Later work |17| . which analyzed dipolar BECs in 
a trap in the limit of lu z /u> p — oo, suggests that a col- 
lapse should still occur even in highly-elongated pancake 
traps, due to a roton-maxon instability involving excita- 
tions with large transverse momenta. As seen from Fig.[7| 
we do observe a collapse of the condensate in the pan- 
cake trap, which, however, is due to the dependence of 
the scattering length on the dipole moment. This col- 
lapse mechanism is similar to that discussed above for 
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FIG. 8: E/N for two particles with hard wall cutoff b = 
0.0137ah o in a spherical trap as a function of the dipole length 
D, obtained by solving the GP equation (solid lines) and 
the Schrodinger equation (dashed lines). The GP line termi- 
nates at points of collapse. For the Schrodinger equation, the 
three lowest eigenvalues are displayed (negative energies are 
not shown). The dotted line shows the results from a GP 
calculation with a constant scattering length a—b. The ver- 
tical dash-dotted lines indicate the resonance positions of V, 
Eq. © (see also Fig.^J. 



spherical traps (see Fig. [3J| . 



D. Excited states and condensate "revival" 

The s- wave-dominated collapse shown in Fig. [31 does 
not necessarily signal the end of the story. After the scat- 
tering length has dropped to -co, it then takes on a large 
positive value (see Fig. EJ) . In this region, the condensate 
is, neglecting three-body recombination effects, perfectly 
stable, at least until the scattering length again becomes 
large and negative. As the resonance is crossed, a new 
two-body bound state appears in the potential Eq. J3J). 
Indeed, this is the first bound state of this potential. 
In the following, we compare the energetics obtained by 
solving the GP equation and the Schrodinger equation for 
N = 2 particles (see Appendix) as the dipole moment is 
tuned across a series of two-body resonances. Extension 
of this analysis to more than N = 2 particles is compli- 
cated by the existence of the two-body bound states of V, 
Eq. (J2J). These two-body bound states imply that solving 
the A-body Schrodinger equation by the DMC method 
requires the use of fixed-node techniques [34(. Such a 
treatment is beyond the scope of this paper. 

Figure [8] shows E/N for two particles under spheri- 
cally symmetric confinement with hardwall cutoff b = 
0.0137a/i O . The exact two-body solution to the 
Schrodinger equation (dashed lines) shows a repeated cy- 
cle of collapse and "revival" . For comparison, solid lines 



8 



show the solution to the GP equation. In the regions 
away from resonance, Fig. [S] indicates good agreement 
between the GP energies and one of the "branches" of the 
exact two-body spectrum. To obtain this agreement, it is 
essential to use the dipole-normalized scattering length, 
which encapsulates the information about the formation 
of bound states of the system, in the effective potential 
that enters the GP equation. Otherwise, if a constant 
scattering length (corresponding to the zero dipole mo- 
ment case with a=b) is used in the GP formalism, the 
unrealistic dotted line in Fig. [S] is obtained. This line 
can only be regarded as a good approximation in certain 
regions between resonances when the dipole dependent 
scattering length is small. As mentioned above, for com- 
putational reasons we have only performed this compar- 
ison for N = 2, but it is strongly suggestive that the 
same should hold for any number of particles. By tuning 
the dipole moment, it should thus be possible to observe 
a collapse followed by a revival of a metastable dipolar 
BEC with increasing dipole strength. The results of a 
GP calculation for a stability diagram of a dipolar BEC 
as a function of the number of particles and the dipole 
moment is shown in Fig. 3 of Ref. [2(j ■ 



In Fig. |21 the resonance positions from Fig. [2] for 
b = 0.0137a/j O are indicated by vertical dash-dotted lines. 
The solution to the two-body Schrodinger equation per- 
sists to slightly larger dipole moments than the reso- 
nance position (this is most apparent for the second 
resonance), because of the added kinetic energy in the 
trap. By contrast, the GP solution collapses at a smaller 
dipole moment, owing to the approximate collapse crite- 
rion (N - l)a(d) = -0.57497a ho . 



At this point we note that our solutions to the N- 
body Schrodinger equation have been restricted to rela- 
tively small number of particles (N < 50). It is some- 
times asserted (see, e.g., Ref. 1]) that the GPE is valid 
only in the limit of large number of particles. This re- 
striction comes from the derivation of the GPE from the 
exact Heisenberg representation equation, by replacing 
the Heisenberg field operator \fr(r) by its classical value 
tya(r). This breaks conservation of particle number, and 
can only be justified for a large number of particles. From 
this point of view, one might be surprised that the mean- 
field equation describes, as demonstrated here, systems 
with small number of particles with high accuracy. How- 
ever, it also possible to derive the GPE from a varia- 
tional ansatz which explicitly conserves particle number 
0,|3j|. This provides the justification for applying the 
GPE even for a small number of particles. In this case, it 
is important to keep the correct N — 1 factor in Eq. JSJ), 
rather than, as is sometimes done, replace it by TV, which 
can only be justified as an approximation valid for a large 
number of particles. 



2.2 
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FIG. 9: Energy per particle E/N versus dipole length D, 
with hardcore radius b — 0.0433a/j O . Solid lines show the GP 
energies for N = 10, 20 and 50 particles (from top to bottom). 
Symbols with error bars show the DMC energies for N = 10 
(dots), N = 20 (circles) and TV = 50 (triangles). 



E. Limits of the mean field approximation 

When the density of the condensate increases, one ex- 
pects the mean field GP theory to become increasingly 
less accurate due to correlations among the particles. In 
the zero-dipole theory of short range interactions the GP 
theory is the lowest order approximation in powers of 
the small parameter na , where n denotes the peak den- 
sity and a the scattering length. The expansion parame- 
ter becomes larger as N and a increase. Treating larger 
number of particles in the DMC simulations increases the 
computational cost. Instead, to check the applicability of 
the GP theory to describe dense dipolar gases, we con- 
sider a larger particle diameter b. Figure shows a com- 
parison between the DMC energies (symbols with error 
bars) and the GP energies (solid lines) for b = 0.0433a^ o 
which is a factor 3.2 larger than we considered in the 
previous sections. A good qualitative agreement still ex- 
ists, but deviations are apparent. The GP calculations 
underestimate the many-body DMC energy, and the de- 
viations grow with increasing number of particles N . The 
GP theory also predicts a collapse at asomewhat smaller 
d than the true value. 

The deviations, which exist even for zero dipole mo- 
ment, are attributed by us mainly to the relatively large 
hard core diameter b = 0.0433a^ o of our "molecules" 
compared to the size of the trap. For zero dipole mo- 
ment, such effects were observed before |36j| . It was found 
that the deviations were reduced by roughly an order of 
magnitude by using a modified GPE that accounts for 
quantum fluctuations [37l l38j . Even better agreement 
was achieved when effective range corrections were taken 
into account [24[. We expect that similar corrections 
would also apply to anisotropic dipolar gases. However, 
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Ref. [36( described condensates with short-range repul- 
sive interactions (a > 0), for which the modified GPE 
included a correction term calculated in the local density 
approximation. In our case there is a long-range dipolar 
interaction which is partly attractive, and furthermore, 
the dipole dependent scattering length becomes negative 
for sufficiently large d. As a result, the local density ap- 
proximation cannot be applied to our problem over the 
whole range, since a homogeneous condensate with at- 
tractive interactions is unstable. 

Finally, we note that additional beyond mean-field cor- 
rections are expected to become important when the 
dipole length £>* becomes larger. Investigation of this 
regime is beyond the scope of this paper. 

IV. CONCLUSIONS 

In conclusion, we have considered a Bose gas with tun- 
able dipolar interactions, controlled by an external field, 
and have shown that taking into account the dependence 
of the scattering length on the dipole moment is essen- 
tial for the correct description of the system within the 
mean field GP method, using the pseudo-potential of Yi 
and You • By comparison with DMC simulations we 
have highlighted both the accuracy of the GP method 
for dipolar gases with low densities, and the existence 
of deviations for dipolar gases with higher densities. The 
mean field theory was shown to predict correctly the con- 
densate size. It also describes well condensates in cylin- 
drical and pancake traps. It was shown that as the dipole 
moment is increased, the scattering length decreases and 
becomes negative as one approaches a resonance corre- 
sponding to the formation of a two-body bound state. 
This effect is the major contribution for the eventual 
collapse of the condensate at a critical dipole moment. 
This is particularly significant for a pancake trap, in 
which case, failing to take the dependence of the scat- 
tering length on the dipole moment into account, col- 
lapse would occur only due to a roton-maxon instabil- 
ity. As the dipole moment is increased further, past the 
two-body resonance, the scattering length a(d) becomes 
positive and we see a "revival" of a stable condensate. 
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APPENDIX 

In this appendix we discuss the numerical techniques 
used in the three different methods of this paper: solution 



to the Gross-Pitaevskii equation, solution to the iV-body 
Schrodinger equation by the DMC method, and solution 
to the two-body Schrodinger equation. 



1. Solution of the Gross-Pitaevskii equation 

Equation JSJ with the pseudopotential Eq. is 
an intcgro-differential equation. The integral term 
d 2 J dr' ^J"^ 9 \ip(r')\ 2 needs special attention due to 
the apparent divergence of the dipolar pseudopotential 
at small distances. The integral does converge if one 
performs the angular integration first. For small dis- 
tances, the density |7/>(r')| 2 may be considered linear 
around r' = r , and integrating the angular part with 
Y20 oc 1 — 3 cos 2 9 gives zero. This can also be seen by 
expanding the density in a multipolc expansion around 
r' = r. Only the Y20 term will contribute to the integral, 
and for a smooth density the coefficient of this term goes, 
for t* — r'| — > 0, as \r — r'\ 2 . Thus, the integral converges 
even without a cutoff. 

Following Ref. 0] , the calculation of the integral can 
be simplified by means of the convolution theorem: 

J dr'V D (r - r')\ip(r')\ 2 = 
F- 1 {F[V](k)FM 2 ](k)} , (A.l) 

where T is the Fourier transform. The Fourier transform 
of the dipolar potential may be performed by expanding 
cxp(zfe -r) in a series of spherical harmonics and spherical 
Bessel functions (the usual expansion of a free planar 
wave in free spherical waves), where only the Y20 term 
gives a non-zero contribution. The result is: 

Ait 

F[V]{k) = -^-(3cos 2 a- I), (A.2) 

where a is the angle between the momentum k and the 
dipole direction. The Fourier transform of Eq. (jA.lfl . 
J^dip] 2 ) is numerically evaluated by means of a standard 
fast Fourier transform (FFT) algorithm, and multiplied 
by J 7 [V](k). Computation of the kinetic energy is also 
accomplished to spectral accuracy through a FFT of the 
wave-function and multiplication by k 2 in momentum 
space, followed by an inverse FFT. 

The ground state of the system was obtained by the 
usual wave-function propagation in imaginary time. The 
propagation was implemented using a 4'th order adaptive 
step Runge-Kutta method. 

This procedure is very computationally intensive. We 
have recently developed an improved method, which we 
also partly used in this work, and which we describe else- 
where. 
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2. Diffusion Monte Carlo simulations 

The many-body Hamiltonian H for a dipolar gas under 
external confinement reads 



H 



N r 

E 



2m J 



1 



j<k 



(A.3) 



where V denotes the interaction potential given by 
Eq. |2"jl. r 3 = (xj,yj,zj) the position vector of the jth 
dipole, r^fe the distance vector, i.e., r^. = — r&, and 
0jfc the angle between the vector and the z-axis. To 
solve the corresponding time-independent Schrodinger 
equation Hip = Eip we employ the variational Monte 
Carlo (VMC) and the Diffusion Monte Carlo (DMC) 
techniques [30j . The former results in a variational bound 
on the energy whereas the latter results in essentially ex- 
act many-body energies. 

To determine variational many-body energies, we write 
the variational wave function ipy as a product of one- 
body terms <p and two-body terms F, 

N N 

i>v{r\,--- ,r N ) = Y[(t>(pj,Zj) x Y[F(r jk ,9 Jk ),(AA) 

j=l j<k 

where <f> contains the Gaussian widths b p and b z , 



4>(p, z) = cxp 



1 / \ 2 1 / \ 2 

1 / p \ 1 I z 



2 V b 



2 V b 



(A.5) 



and F the hardcore radius b and the "shape" parameters 
pi through p 5 , 



F(r,i 



1 



1 + ^7 (pi+P2 cos 



Pi 

r P5 



cos 4 e 



;a.6) 



To first order the parameters of p are determined by 
the external confining potential and those of F by the 
two-body potential V . For a non-interacting gas with 
b = d = 0, e.g., the variational wave function ipy with 
b p = b z = 1 and p\ = P2 = Pa = is an exact solution to 
the many-body Schrodinger equation. For a finite hard- 
core radius 6, the term in the first pair of square brackets 
on the right hand side of Eq. I)A.6(I coincides with the 
low-energy s-wave scattering wave function for two par- 
ticles interacting through a spherically symmetric hard- 
core potential with radius b, i.e., this term ensures that 
the many-body wave function goes to zero as the dis- 
tance r between two particles approaches the hardcore 
radius b. The functional form of the angle-dependent 
term of F [term in the second pair of square brackets 
on the right hand side of Eq. (|A.6|l ] is motivated by the 
shape of the essentially exact two-body wave functions 



for non- vanishing d 2 , which is obtained numerically us- 
ing B-splines (see Appendix, subsection l"3*)l . 

The variational parameters b p , b z and p\ through p$, 
collectively denoted by p, are optimized for each N, u> p , 
lo z and -D* by minimizing the energy Ey, 



E V (P) = 



^v(p)\H\^y(p)) 
(ipv(pMv(p)) 



(A.7) 



where the integration over the 3N coordinates r%, ■ ■ ■ , rjy 
is performed using Metropolis sampling. Our 
parametrization of the variational wave function ipy pro- 
vides an excellent description of weakly-interacting gases 
with small nb 3 and nD%, where n denotes the condensate 
density at the trap center. As b/aho or D^/aho increase, 
i.e., as the gas becomes more strongly- interacting, the 
VMC energy Ey recovers an increasingly smaller frac- 
tion of the essentially exact DMC energy. 

The optimized variational wave function ipy is sub- 
sequently used as a guiding function in our DMC cal- 
culations with importance sampling. The DMC tech- 
nique solves the many-body Schrodinger equation for 
the ground state energy and for structural properties 
by starting with an initial "walker distribution" , which 
can be thought of as a stochastic representation of the 
many-body wave function, and by then projecting out 
the lowest stationary eigenstate through propagation in 
imaginary time. To treat high-dimensional systems (here 
with as many as 150 degrees of freedom), the short-time 
Green's function propagator is evaluated stochastically. 
The resulting DMC energies are essentially exact, i.e., 
they are independent of ipy, with the only uncertainty 
stemming from the finite time step used in the propa- 
gation. For the DMC energies reported, the time step 
errors are smaller than the statistical uncertainties. Us- 
age of a good guiding function, i.e., our optimized ipy, is 
essential for the DMC algorithm to be numerically sta- 
ble. If ipy coincides with the exact many-body wave 
function — as is the case for non-interacting inhomoge- 
neous gases — , the resulting DMC energies have vanish- 
ing variance and thus vanishing errorbars. Interactions, 
and correspondingly approximate ipy, introduce statis- 
tical uncertainties in the DMC energies, which can be 
reduced by increasing the computational efforts. We cal- 
culate structural expectation values using a descendant 
weighting scheme [39l l40l El) , which in principle elimi- 
nates any dependence of the structural expectation val- 
ues on ipy. However, our structural expectation values 
may be slightly biased by ipy since there is a tradeoff be- 
tween the statistical uncertainty and the lengths of the 
"side walks" used in the descendant weighting scheme. 

The standard DMC algorithm outlined above describes 
the ground state. To describe metastable condensates, 
i.e., excited many-body states with gas-like character, our 
DMC simulations take advantage of the topology of the 
underlying configuration space. The metastable conden- 
sate state, characterized by large interparticle distances, 
is separated from bound many-body states, character- 
ized by small interparticle distances, by a "barrier" (see 
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Section Till Af l . In this respect, the topology of the high- 
dimensional configuration space at hand is similar to the 
configuration space of a one-dimensional double-well po- 
tential with large barrier. For large enough barrier, the 
two regions in configuration space correspond to two ef- 
fectively orthogonal Hilbert spaces. This implies that 
DMC walkers with initial coordinates corresponding to a 
metastable gas-like state have a vanishingly small proba- 
bility to tunnel into the region of configuration space cor- 
responding to molecular-like bound states and that the 
simulations consequently converge to the metastable con- 
densate state and not to energetically lower-lying cluster 
states. The presence of the barrier is thus crucial for de- 
scribing metastable condensate states with negative scat- 
tering length by the DMC algorithm. 



3. Numerical solution of Schrodinger equation for 
two dipoles 



For two dipoles in a trap interacting through the po- 
tential of Eq. (J2J) , the Schrodinger equation is separable 
in relative distance and center of mass coordinates. The 
center of mass equation describes the harmonic motion 
of the two dipoles as a whole in the trap, and has the 
usual harmonic oscillator ground state. The relative dis- 
tance equation has cylindrical symmetry and was solved 
numerically by expanding the wavefunction on a basis 
set of two-dimensional B-splines, with the appropriate 
boundary conditions. The full Hamiltonian matrix in 
this basis was constructed and diagonalized. 
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